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INTRODUCTION 


A  substantial  need  has  been  demonstrated  within  Armament  Research  and 
Development  Center  (ARDC)  for  a  convenient— to— use  and  computationally  efficient 
flight  simulation  for  smart  munitions.  The  simulation  must  be  easily  adaptable 
to  varying  projectile  designs  and  especially  to  changes  in  autopilots  which 
undergo  rapid  evolution.  During  the  design  cycle  of  the  projectile,  a  large 
number  of  individual  flight  simulations  must  be  performed.  As  an  example,  during 
the  development  of  the  Copperhead  projectile,  the  ARDC  scientific  computers  were 
overwhelmed  by  the  workload  and  a  substantial  amount  of  the  available  time  on  the 
BRL  CDC  7600  supercomputer  was  dedicated  to  the  computations.  In  spite  of  the 
continuous  development  and  improvement  in  computer  systems,  flight  simulations 

will  continue  to  provide  a  substantial  workload  which  must  be  accomplished  effi- 
ciently . 

To  meet  this  need,  ARDC  is  formulating  methods  and  writing  modular  software 
for  rapidly  developing  autopilot  simulations  for  guided  projectiles  and  muni¬ 
tions.  An  innovative  technique  was  developed  by  which  exact  analytical  solutions 
to  the  required  transfer  functions  are  applied  in  a  piecewise  manner  within  a 
larger  but  lower  frequency  problem  which  must  be  solved  numerically. 

Since  the  time  constants  associated  with  the  autopilot  components  are  often 
small  compared  with  their  driving  terms,  the  integration  time  step  is  driven  to 
very  small  values  to  obtain  stable  numerical  integrations  which  results  in  very 
long  computer  run  times.  The  use  of  piecewise  analytical  solutions  to  the  trans¬ 
fer  functions  guarantees  valid  integration  of  the  autopilot  transfer  functions 
regardless  of  the  integration  time  step,  provided  that  the  driving  terms  vary 
inappreciably  over  the  time  step.  This  is  a  less  stringent  requirement  than  is 
needed  for  stable  numerical  integration  since  the  driving  term  rates  are  commen¬ 
surate  with  the  airframe  time  constants,  which  are  typically  large  compared  to 
that  of  the  autopilot. 

Five  transfer  functions  were  solved  in  closed  form  to  provide  a  fast 
executing  computer  code.  These  particular  transfer  functions  were  selected  since 
an  autopilot  can  generally  be  represented  by  a  concatenation  of  these  functions 
to  rate  sensors,  switches,  limiters,  and  dead  zones. 


DISCUSSION 


Apart  from  switches,  limiters,  and  dead  zones,  autopilot  transfer  functions 
are  Laplace  transform  representations  of  differential  equations.  Incorporated 
into  a  digital  simulation,  transfer  functions  are  converted  into  differential 
equations  and  are  usually  solved  numerically.  Sometimes  this  causes  considerably 
longer  run  times  than  simulation  of  unguided  projectiles. 

The  time  step  required  to  perform  the  integration  is  driven  by  two  con- 
Th^  driving  term  or  input  must  not  vary  appreciably  in  the  time 
step.  (2)  The  time  step  must  be  sufficiently  small  to  Insure  a  stable 
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integration.  If  the  time  constants  associated  with  the  autopilot  are  small  com¬ 
pared  to  those  of  the  airframe  (which  is  typical),  the  integration  time  step  will 
be  driven  to  very  small  values  and  run  times  will  be  very  long.  Since  the 
driving  term  is  a  product  of  the  airframe  motion,  an  inherently  slower  process, 
stable  integration  is  a  stronger  driver  to  fine  integration  than  the  driving 
term's  remaining  essentially  constant  during  the  Integration  time  step.  By 
analytically  solving  these  transfer  functions  stepwise  with  constant  driving 
terms,  the  second  constraint  is  eliminated.  The  size  of  the  integration  time 
step  is  limited  only  by  the  first  constraint.  Therefore,  the  use  of  closed  form, 
analytic,  stepwise  solutions  subject  only  to  the  first  constraint  can  lead  to 
larger  integration  time  steps  and  shorter  run  times. ^ 

Only  five  transfer  functions  are  needed  to  handle  the  typical  autopilot: 
the  first-order  lag,  the  first-order  lag  with  differentiation,  the  fir-  -order 
lead/lag,  the  first-order  lag  with  integrator,  and  the  second-or  .er  lag/ 
oscillator.  (Guidance  and  control  systems  often  contain  rate  sensors.  The 
treatment  of  rate  sensors  will  be  published  in  a  separate  report.)  These 
transfer  functions  are  described  and  solved  analytically  below;  the  computer  code 
implementing  the  solutions  is  shown  in  appendix  A.^ 

"^^^^^Shout  this  development,  the  driving  term  is  assumed  not  to  vary  appre¬ 
ciably  during  a  time  step.  Care  must  be  taken  to  adjust  the  time  step  downward 
as  necessary  to  satisfy  this  requirement.  To  take  advantage  of  this  technique 
for  reducing  run  time,  the  time  step  should  be  increased  whenever  the  driving 
terras  are  varying  slowly. 


ANALYSIS 


First-Order  Lag 


•A  first-order  lag  is  represented  by  the  Laplace  operator  s  as 
l/(Ts  +  1) 

In  differential  equation  form 
T  dy/ dt  +  y  =  D 


Michael  J.  Amoruso  and  Dennis  D.  Ladd,  "TELUM,  A  Comprehensive  Digital  Flight 
Simulation  of  the  Copperhead  Projectile,"  Special  Publication  ARLCD-SP-82003 , 
U.S.  Army  Armament  Research  and  Development  Command,  Dover,  NJ,  June  1982. 

Analytical  methods  used  in  this  study  were  derived  from  John  J.  D'Azzo  and 
Constantine  H.  Houpis,  Feedback  Control  Systems  Analysis  and  Synthesis, 
McGraw-Hill,  New  York,  1960;  and  Murray  R.  Spiegel,  implied  Differential 
Equations .  Prentice-Hall,  Englewood  Cliffs,  NJ,  1958. 
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where 


y  -  output  or  dependent  variable 
t  =  time  or  independent  variable 
T  =  time  constant 


D  =  driving  term 

The  general  solution  to  equation  1  is 
general  solution  to  the  homogeneous  form, 
by  setting  D  to  zero  and  can  be  written  as 


the  sura  of  a  particular  solution  and  a 
The  homogeneous  equation  was  obtained 


-  T  dy/y  =  dt 


which  can  be  verified  to  yield  the  solution 


where^A  is  the  Integration  constant  determined  from  the  initial  conditions  y  =  y 
geneo^s  ‘iq’uatl' n  conditions  yields  the  general  solution  to  the  ho«o5 


-(t  -  t  )/T 

y  =  y^e 


(3) 


tion^  particular  solution  to  equation  1,  the  expression  is  verified  by  substitu- 


y  = 


1 


(4) 


Adding  equations  3  and  4  gives  the  complete  general  solution 


y(t) 


+  1 


(5) 


Note  that  equation  5  has  the  correct  limit  of  D  as  T  goes  to  zero 
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^rst— Order  Lag  with  Differentiation 


This  transfer  function  has  the  Laplace  operator  representation 


s/(Ts  +  1) 


or  the  time  domain  form 


T  dy/dt  +  y  =  dD/dt 


(6) 


This  representation  is  not  desirable  since  the  driving  term  appears  as  a  time 
derivative  that  would  generally  have  to  be  evaluated  numerically.  However,  equa¬ 
tion  6  can  be  recast  as  a  pair  of  coupled  differential  equations  that  do  not 
contain  any  derivatives  of  the  driving  term  D.  An  auxiliary  variable  z  is  intro¬ 
duced  and  the  equation  is  replaced  by  the  expressions 


dz/dt  =  y 


(7) 


and 


y  =  (D  -  z)/T 


(8) 


Differentiating  equation  8  and  using  equation  7  to  eliminate  dz/dt  verifies  that 
these  two  equations  are  equivalent  to  equation  6.  Note  that  the  derivative  of  D 

does  not  appear.  By  combining  equations  7  and  8,  the  expression  for  z(t)  is 
obtained. 


dz/dt  +  z/T  =  D/T 


(9) 


and,  after  integration,  the  results  for  z(t)  substituted  into  equation  8.  Equa¬ 
tion  9  can  be  integrated  from  equation  5  by  noting  the  identity  of  equations  1 
and  9.  Therefore, 


-(t  -  t J/T 

y(t)  =  e  (d  -  Z^)/T 


(10) 


where  =  2  at  t 
equation  8. 


to*  The  auxiliary  variable  z^  can  be  eliminated  by  using 


h 


y(t) 


\ 


(d  -  D  )  +  y 
O  o 


-(t 


(11) 


Note  that  the  exponential  over  T  goes  to  l/(t  -  t  ) 
equation  11  has  the  correct  limit  dD/dt  as  T  goes  to 


by  L' Hospital's  rule  so  that 
zero. 


First-Order  Lead/Lag 


This  transfer  function  is  represented  by 


[T2S  +  1)/[TiS  +1)  0 


or 


Tj  dy/dt  +  y  =  T2  dD/dt  +  D 


(12) 


The  derivative 
is  eliminated 
equations 


w  f'l  driving  term  is  usually  not  available  in  analytic  form.  It 
by  introducin::  an  auxiliary  variable,  z,  and  a  pair  of  coupled 


dz/dt  =  (D  -  y)/Tj 


(13) 


and 


y  =  [T2/TJ  D  +  z 


(14) 


The  equivalence  with  equation  14 
^fimlnating  dz/ dt  with  equation 
equation  14  as  follows: 


can  be  verified  by  differentiating  it  and  then 
13.  To  obtain  an  expression  for  y(t),  integrate 


T^  dz/dt  +  2  =  (1  -  T2/TJ  D 


(15) 


5 


This  is  obtained  by  substituting  equation  14  into  equation  13,  and  then  substi¬ 
tuting  the  result  for  z(t)  into  equation  14  to  obtain  y(t).  The  explicit  solu- 
not  be  written*  Instead,  an  algorithm  will  be  provided  based  on  the 
observation  that  equation  15  is  equivalent  to  equation  1  if  the  following  substi¬ 
tutions  are  made: 


Ti  T  z  y  (l  -  T2/TJD  D  (16) 

Algorithm: 

1.  Substitute  equation  16  into  equation  5  to  obtain  z(t). 

2.  Save  z(t)  for  use  as  z^  in  the  next  integration  step. 

3.  Substitute  z(t)  into  equation  14  to  obtain  y(t). 

Note  that  the  lag  time  constant  cannot  be  zero  in  equation  16,  but  a  vanish¬ 
ing  lead  time  constant  yields  the  correct  limiting  case  of  a  simple  first-order 
lag. 

First-Order  Lag  with  Integrator 

The  representation  of  this  transfer  function  is 
1/s  (Ts  +  1) 


or 


T  d^y/dt^  +  dy/dt  =  D 


(17) 


Substitution  verifies  that  the  following  expression  is  a  particular  solution  to 
the  nonhomogeneous  differential  equation  17  since,  by  assumption,  D  is  constant 
during  the  Integration  time  step: 


y  =  Dt 


(18) 


The  homogeneous  form  of  equation  17  can  be  written  as  the  coupled  equations 


T  dz/dt  +  z  =  0 


(19) 


and 


z  =  dy/dt 


(20) 


Equation  19  is  similar  to  equation  1  with  the  substitution  of  z  for  y  and 
or  D;  therefore,  the  solution  can  be  obtained  from  equations  2  and  20 


zero 


dy/dt  =  z(t)  = 


or 


y(t) 


=  -Te"^*"  ■ 


+  B 


(21) 


tlL7r^^9n  to  equation  17  is  the  sum  of  the  general  homogeneous  solu¬ 

tion  (,eq  21)  and  the  particular  nonhomogeneous  solution  (eq  18) 


y(t) 


+  Dt  +  B 


(22) 


dy/dt  =  e  +  D  (23) 

The  integration  constants  A  and  B  can  be  determined  by  invoking  the  initial  con¬ 
ditions  y  =  y^  and  dy/dt  =  y'  at  t  =  t^ 


y 


o 


a)/t 

+  D  t  +  B 
o 


(24) 


y’ 


-(t„  -  a)/t 

e  +  D 


(25) 


The  integration  constant  A  can  be  eliminated  from  equation  24 
25  to  eliminate  the  exponential  term  and  solving  for  B 


by  using  equation 
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(26) 


B  =  Yo  +  T  (y'  -  D)  -  Dt^ 

Taking  the  natural  logarithm  of  equation  25  yields 
A  =  T  In  (y’  -  D)  + 

The  explicit  solution  is  therefore 


and 


y(t)  =  y^  +  T  (y'  -  D) 


+  D  ft  -  t  ) 


(27) 


(28) 


-(t  -  t  )/T 

dy/dt  =  (y'  -  D)  e  °  +  D 


(29) 


Note  that  equations  28  and  29  go  to  the  correct  limits  if  T  goes  to  zero 


Second-Order  Lag/Oscillator 

This  transfer  function  can  be  represented  in  the  form 


l/(Is2  +  Ds  +  K)  K  *  0,  1  *  0 

or 

I  d^y/dt^  +  dy/dt  +  Ky  =  T 

The  homogeneous  solution  to  equation  30  can  be  verified  to  be 
y(t)  =  e^^ 

by  substituting  into  the  homogeneous  form  of  equation  30  to  obtain 
istic  equation 

1L2  +  dL  +  K  =  0 


(30) 


(31) 

the  character— 

(32) 
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which  has  the  following  roots : 


r 

M  =  !  -D  + 


(33) 


(34) 


The  following  three  cases  are  treated  separately  according  to  whether  the  radi- 
cand  IS  positive,  negative,  or  zero. 


Case  1.  Positive  radicand,  damped  solution 

The  homogeneous  and  nonhomogeneous  particular  solutions  can  be  verified  to 


yg(t)  =  Ae'^*'  +  Be^^ 


(35) 


yp(t)  =  T/K 


(36) 


The  complete  solution  is 


y(t)  =  Ae^*^  +  Be’^^  +  T/K 


(37) 


dy/dt  =  MAe^^  +  NBe^*" 


(38) 


The  initial  conditions  y 
A  and  B.  ^ 


y  at  t^  and  y'  -  dy/dt  at  t^  determine  the  constants 


Mt  Nt 

y^  =  Ae  °  +  Be  °  +  T/K 


(39) 


9 


Mt  Nt 

y'  =  MAe  +  NBe  ° 

Solving  simultaneously  gives 


-N  t 


=  (y^  -  y’/M  -  T/K)e  °/  (l  -  N/M) 


(40) 


(41) 


-M  t 

A  *  (y^  -  y'/N  -  T/K)e  °/  (l  -  M/n) 


(42) 


Case  2.  Zero  radlcand,  critically  damped 

The  characteristic  equation  has  only  one  distinct  root 


L  =  -D/21 


(43) 


The  homogeneous  general  solution  and  nonhomogeneous  particular 
verified  to  be 


solution  can  be 


/*.\  .  Lt  Lt 

yjj  (t)  =  Ae  +  Bte 


(44) 


yjj(t)  -  T/K 


(45) 


The  complete  solution  Is 


y(t)  =  (A  +  Bt)  e^^  +  T/K 


(46) 


dy/dt  =  (LA  +  B  +  LBt)e^*^ 


(47) 
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Using  the  Initial  conditions  y  »»  y  at  t  and  y' 
constant  can  be  determined  as  follows:  ° 


dy/dt  at  t^,  the  integration 


Lt 

=  (a  +  Bt^)e  °  +  T/K 


(48) 


Lt 


y'  =  (la  +  B  +  LBt  )e  ° 


(49) 


-Lt 


B  = 


=  e  [y'  +  L  (t/K  -  y^)] 


(50) 


-Lt 


A  =  [(y^  -  t/K)  (1  +  L  t^)  -  t^  y']e  ° 


(51) 


Case  3.  Negative  radicand,  oscillation 

The  roots  of  the  characteristic  equation  become  complex  numbers 


0)  = 


(K/I)  -  (D/21) 


V; 


(52) 


The  homog:;neous  and  nonhomogeneous  particular  solutions 


are 


Lt 


yjj(t)  =  Ae  sin  (tot  +  ()>)+  T/K 


(53) 


yH(t)  =  T/K 


(54) 
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where  A  and  ({>  are  integrati 


on  constants.  The  complete  solution  is 


y 


Ae 


Lt 


sin  (ojt  +  (())+  T/K 


(55) 


dy/dt  -  LAe  sin  (wt  +  (j))  +  wAe^^  cos  (tot  +  ({>) 


(56) 


Using  the  initial 
gives 


conditions  and  y’  and  solving  for  A  in  the  expression  ror  y 


-L  t 

^  ~  [Yq  “  T/K]e  °/sin  (tot^  +  <))) 


(57) 


where 


sin  (tot^  +  (())  0 


(58) 


If  equation  58  fails,  A  can  be  obtained  from  the  expression  for  y' 


-L  t 

A  -  y'e  /[l  sin  [u  t^  +  ij))  +  q)  cos  +  (})]] 


(59) 


Since  0)  in  equation  52  is  always  positive  definite  and  the  sine  and  cosine  are 

slo^s^Ttn’l  ^  59  is  always  defined.  Note  that  the  expres- 

whethpr  '*’*  Solving  for  <j)  yields  the  same  result 

57  is  substituted  into  the  Initial  condition  on  dy/dt  or  equa¬ 
tion  59  into  the  initial  condition  on  y. 


<(,  = 


-1 

tan 


0)  (y^  -  t/k] 

[y'  -  i'(yo  ~  t/k)] 


(0  t 

o 


(60) 


E^LtlJn^fiff  in  equation  57  or  59  by  substituting  equation  60. 

aquation  60  shows  that  equation  58  is  equivalent  to 


y^  -  T/K  =  0 


(61) 
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The  imposition  in  equation  30  of  the  requirement  that  parameters  K  and  I  not 

3^  45  Tnd^'sT^r'''  ^^lidity  of  the  solutions.  For  example,  equations 
36,  45,  and  54  become  singular  when  K  vanishes.  Similarly,  the  roots  of  the 
characteristic  equation  are  singular  if  I  is  zero  (eqs  33,  45,  and  52).  This  is 
not  a  limitation  since  equation  30  reduces  to  a  first-order  lag  with^ integrator 
when  K  vanishes  and  to  a  simple  first-order  lag  when  I  vanishes. 


However,  the  parameter  D  in  equation  30  may  be  zero, 
characteristic  equation  32  degenerates  to 


If  D  vanishes,  the 


IL2  +  K  =  0 


M  = 


±  i  0) 


Therefore,  the  roots  are  imaginary  corresponding  to  a  pure  oscillation  with  no 
amping.  Only  case  3  should  be  considered.  Nowhere  does  D  appear  in  the  denomi¬ 
nator  of  any  expression;  therefore,  the  solution  remains  valid  if  D  is  zero. 


RESULTS 


on  thl^%n^  saving  resulting  from  this  piecewise  analytical  technique  depends 
^  execution  of  the  competing  numerical  method.  A  sample  run  is 

ncluded  in  appendix  B  so  that  comparisons  can  be  made  with  other  preferred 
numerical  integration  techniques.  For  comparison,  the  last  case  (KTFId'iI)  was 
run  analytically  with  time  steps  of  0.01  and  0.005  second  and  run  numerically 
using  Mvanced  Continuous  Simulation  Language  (ACSL).^  whereas  the  analytical 
method  gave  the  same  results  for  time  steps  0.001  and  0.005  second,  ACSL  could 
not  operate  with  a  time  step  of  0.001  second.  ACSL  was  able  to  obtain  nearly 

second^  The^se  r'  ^  the  third 

second.  These  results  are  compared  as  follows: 


Time 

(sec) 


_ Analytical _ 

DT=0.01  sec  DT=0.005  sec 


Numerical  (ACSL) 
DT=0 .00001  sec 


0  10.00000 

1  1.54640 

2  -9.52188 

3  -4.49206 


10.00000 

1.54640 

-9.52188 

-4.49206 


10.00000 

1.54381 

-9.52327 

-4.48494 


^vanced  Continuous  Simulation  Language,  User  Guide/Reference  Manual.  MitPh.n 
and  Gauthier  Associates,  Inc.,  Concord,  MA,  1981.  - - - 
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In  this  particular  case,  the  analytical  technique  resulted  in  at  least  two  orders 
of  magnitude  improvement  in  integration  step  size. 


The 

result  is 


GONCLUSIONS 

Y 


modular,  piecewise  analytical  technique  presented  in  this  report  can 
considerable  savings  in  analysis  and  associated  computer  time. 
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APPENDIX  A 

PROGRAM  LISTING  AND  SAMPLE  EXECUTION 
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o  o 


PROGRAM  LISTING 


PROGRAM  TFMAIN ( INPUT , OUTPUT , TAPE5= INPUT , TAPE6=OUTPUT ) 

MAIN  PROGRAM  TO  DRIVE  TRANSFER  FUNCTION  SUBROUTINE. 

FOR  TEST  PURPOSES  ONLY. 

COMMON/TRNSF/TAUl (50)  ,TAU2(50)  ,TAU3(50)  ,TFTIME(50)  ,TFAUX  (50 )  , 
+  TFAUX2  (50)  ,KTFTYP (50)  ,GAIN (50) 

OPEN  (8) 

WRITE  (8,11) 

C  READ  INPUT  DATA  AND  INITIALIZE. 

CALL  TRNF  10 
WRITE  (8,11) 

C  200  DEFINE  INTEGRATION  TIME  STEP  DT . 

DT  =  .005 
DT2=  DT/2. 

WRITE(8,13)  DT 
DO  2  J=l,13 
WRITE  (8,12) 

C  100  DEFINE  INITIAL  TIME  T  AND  FINAL  TIME  TMAX . 

T  =  0.0 

TMAX  =  5.0 

KMAX  =  IFIX (TMAX/DT) +1 
DO  1  K=0,KMAX 
T  =  DT*FLOAT(K) 

C  300  DEFINE  DRIVING  TERM  "DRIVE" 

,  DRIVE  =  COS (2.*T)  -  i.O 

CALL  TRNF  (DRI VE , XOUT , T , J ) 

IF (ABS (IFIX (T) -T) .LT.DT2)  WRITE(8,10)  J , T , DRI VE , XOUT 

1  CONTINUE 

2  CONTINUE 
STOP 

10  FORMAT(1X,I5,3F20.5) 

11  FORMAT (’ 1 ' ) 

12  FORMAT(//,'  LTFID' ,16X, 'TIME'  ,14X, 'DRIVE' ,19X,  'X' ) 

13  FORI'.AT  (  '  THE  INTEGRATION  TIME  STEP  DT  =  ',F9.6,'  '  //) 

’^ND  * 


C 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  TRNF  10 

LTFTYP  1  ->  GAIN/(TAU1*S+1) 

2  ->  GAIN*S/(TAU1*S+1) 

3  ->  GAIN/(S* (TAU1*S+1) ) 

4  ->  GAIN* (TAU2*S+1)/(TAU1*S+1) 

5  ->  GAIN/(TAU1*S**2  +  TAU2*S  +TAU3) 

LTFID  ID  NO.  OF  CHANNEL  (VARIABLE) . 

TUI  FIRST  TIME  CONSTANT 

TU2  SECOND  TIME  CONSTANT  (LTFTYP=4  OR  5) 

TU3  THIRD  TIME  CONSTANT  (LTFTYP=5) 
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o  o  n 


XO 

INITIAL 

CONDITION 

ON 

variable 

DXO 

INITIAL 

CONDITION 

ON 

DERIVATIVE  (LTFTYP=5) 

AA 

ALPHAMERIC  LABEL 

FOR 

PRINTOUT 

COMMON/TRNSF/TAUl (50)  ,TAU2 (50)  ,TAU3  (50)  ,TFTIME (50)  ,TFAUX  (50)  , 
+  TFAUX2  (50)  ,KTFTYP (50)  ,GAIN (50) 

CHARACTER  AA* 30 , LABEL ( 5 ) *30 
DATA  LABEL/'GAIN/(TAU1*S+1) 

+  'GAIN*S/(TAU1*S+1) 

+  'GAIN/(S* (TAU1*S+1) ) 

+  'GAIN* (TAU2*S+1)/(TAU1*S+1) 

+  'GAIN/(TAU1*S**2+TAU2*S+TAU3)  '/ 


WRITE  (8,1050) 

DO  1  1=1,50 

TAUl(I) 

0. 

TAU2(I) 

0. 

TAU3(I) 

0. 

TFTIME(I) 

0. 

TFAUX(I) 

0. 

TFAUX2(I) 

0. 

KTFTYP(I) 

0 

GAIN(I) 

0. 

1  CONTINUE 
IFLAG  =  0 
DO  2  1=1,50 

READ  (5,1051,END=4)  LTF ID , LTFTYP , TUI , TU2 , TU3 , GA , XO , DXO , AA 
IF  (LTFTYP.lt. 1  .OR.  LTFTYP. GT. 5)  THEN 
WRITE (8 ,1052) LTFID,LTFTYP,TU1,TU2,TU3,GA,X0,DX0 ,AA 
IFLAG  =  1 
WRITE (8,1004) 

GO  TO  2 
ELSE 

WRITE (8,1052)LTFID,LTFTYP,TU1,TU2,TU3,GA,X0,DX0,AA,LABEL (LTFTYP) 
END  IF 

IF  (LTFID.lt. 1  .OR.  LTFID.GT.50)  THEN 
IFLAG  =  1 
WRITE  (8 , 1001) 

GO  TO  2 
END  IF 

IF (TAUl (LTF ID)  .NE. 0. .OR.TAU2 (LTF ID)  .NE. 0. .OR.TAU3  (LTF ID)  . NE . 0 . 

+  .OR.KTFTYP  (LTFID)  .NE.O)  THEN 
IFLAG  =1 
WRITE (8,1002) 

GO  TO  2 
END  IF 

IF(TU1.GT.0. .AND.TU2.EQ.0. . AND . TU3 . GT . 0 . . AND . LTFTYP . EQ . 5 ) 

+  GO  TO  3 

IF  ( (TUl.LE.O. .OR.TU2.LE.O. .OR.TU3.LE.O.) . AND . LTFTYP . EQ . 5 ) THEN 
WRITE (8 , 1003 ) 

IF  (TUl.EQ.O.  .AND.  TU3.GT.0.  . AND . TU2 . GT . 0 . )  THEN 

WRITE (8 , 1005) 

LTFTYP  =  1 
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GA  =  GA/TU3 
TU1=TU2/TU3 
TU2=0 . 

TU3=0 . 

WRITE  (8,1000) 

WRITE (8,1052) LTFID,LTFTYP,TU1,TU2,TU3,GA,X0,DX0,AA 
+  , LABEL (LTFTYP) 

GOTO  3 
END  IF 

IF  (TU3.EQ.0.  .AND.  TUl.GT.O.  . AND . TU2 . GT . 0 . )  THEN 
WRITE (8 , 1006) 

LTFTYP  =  3 
GA  =  GA/TU2 
TU1=TU1/TU2 
TU2=0 . 

TU3=0. 

WRITE (8 , 1000) 

WRITE  (8 ,1052)LTFID,LTFTYP,TU1,TU2,TU3,GA,X0,DX0,AA 
+  , LABEL (LTFTYP) 

GOTO  3 
END  IF 

IFLAG  =  1 
GOTO  2 
END  IF 

IF  (TUl.LE.O.  .AND.  LTFTYP. LE. 4)  THEN 
IFLAG  =1 
WRITE (8,1003) 

GO  TO  2 
END  IF 
3  CONTINUE 

IF  (TU2.GT.0.  .AND.  LTFTYP. LE. 3)  THEN 
TU2=0 . 

WRITE (8,1003) 

END  IF 

IF  (LTFTYP. EQ. 4  .AND.  TU2.EQ.0.  .AND.  TUl.GT.O.)  THEN 
LTFTYP  =  1 
WRITE  '8 , 1003 ) 

WRITE (8 , 1005) 

WRITE (8,1000) 

WRITE (8 , 1052) LTF ID, LTFTYP, TUI, TU2,TU3 ,GA,X0 ,DX0 ,AA 
+  , LABEL (LTFTYP) 

END  IF 

IF  (TU1.EQ.TU2  .AND.  LTFTYP. EQ. 4)  THEN 
WRITE (8,1007) 

END  IF 

IF  (TU3.GT.0.  .AND.  LTFTYP. LT. 5)  THEN 
TU3=0 . 0 
WRITE (8,1003) 

END  IF 

IF  (GA.LE.O.)  WRITE (8,1008) 

TAUl  (LTFID)  =  TUI 
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1000 

1001 

1002 

1003 

1004 

1005 

1006 


FUNCTION  WILL  BE  TREATED  AS: ' 
10***  KTFID  OUT  OF  RANGE.'  ) 
10***  THIS  TRANSFER  FUNCTION 


FORMAT ( ' 
FORMAT  ( ' 
+  'KTFTYP 
1007  FORMAT (' 


TAU  OUT  OF  RANGE.'  ) 
KTFTYP  OUT  OF  RANGE.') 
LAG.  KTFTYP  ->  1.' 
LAG  WITH  INTEGRATOR. 

TRANSFER 


RANGE. 


) 


KTFID 


'GAIN' ,9X, 'X 


TAU2(LTFID)  =  TU2 
TAU3 (LTFID)  =  TU3 
TFAUX(LTFID)  =  XO 
TFAUX2 (LTFID) =  DXO 
KTFTYP  (LTFID)  =  LTFTYP 
GAIN (LTFID)  =  GA 
WRITE (8,1053) 

WRITE (8,1053) 

IF  (IFLAG.NE.O)  STOP 
RETURN 

FORMAT ( '  ***THIS  TRANSFER 
FORMAT ('  ***ERROR  IN  TRNF 

FORMAT ('  ***ERROR  IN  TRNF 

+  , 'ALREADY  DEFINED.') 

FORMAT ('  ***ERROR  IN  TRNF  10*** 

FORMAT ('  ***ERROR  IN  TRNF  10*** 

SUBSTITUTING  1ST  ORDER 
SUBSTITUTING  1ST  ORDER 
->  3.  '  ) 

***ERROR  IN  TRNF  10***  TAUl  =  TAU2 
+  '  DEGENRATES  TO  UNITY.'  ) 

1008  FORMAT ('  ***ERROR  IN  TRNF  10***  GAIN  OUT  OF 

1050  FORMAT (///,' ITRANSFER  FUNCTIONS  DEFINITIONS 
+  '  KTFTYP' ,6X, 'TAUl' ,6X, 'TAU2' ,6X, 'TAU3' ,6X 
+  8X, 'DX' ,//) 

1051  FORMAT  (212, 5X,6F6.0,5X,A30) 

1052  FORMAT (217, 6F10. 5,/, '  DESCRIPTION:  ' , A30 , 2X , A30 ) 

1053  FORMAT (/) 

END 

SUBROUTINE  TRNF (XIN , XOUT , T , KTF ID ) 

C  THIS  ROUTINE  EVALUATES  THE  FOLLOWING  TRANSFER  FUNCTIONS 

C  1ST  ORDER  LAG; 

C  1ST  ORDER  LAG  WITH  DIFFERENTIATION; 

C  1ST  ORDER  LAG  WITH  INTEGRATION; 

C  COMBINED  LEAD/LAG; 

C  2ND  ORDER  LAG  -  HARMONIC  OSCILLATOR. 

C  XIN  DRIVING  TERM. 

C  XOUT  OUTPUT  OF  TRANSFER  FUNCTION. 

C  T  TIME. 

C  KTFTYP  1  ->  1/(TAU1*S+1) 

C  2  ->  S/(TAU1*S+1) 

C  3  ->  1/(S* (TAU1*S+1) ) 

C  4  ->  (TAU2*S+1) / (TAU1*S+1) 

C  5  ->  1/(TAU1*S**2  +  TAU2*S  +TAU3) 

C  KTFID  ID  NO.  OF  CHANNEL  (VARIABLE) . 

COMMON/TRNSF/TAUl (50)  ,TAU2 (50)  ,TAU3  (50)  ,TFTIME (50)  ,TFAUX (50)  , 
+  TFAUX2  (50)  , KTFTYP  (50)  , GAIN (50) 

REAL  KK,II 

IF  (KTFID.lt. 1  . OR. KTFID. GT. 50)  GOTO  102 
KTF  =  KTFTYP (KTFID) 

IF  (KTF.LT.l  .OR.  KTF.GT.5)  GOTO  103 
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XOUT  =  0. 

XTMP  =  0. 

XTMP2  =  0. 

IF  (TFTIME (KTFID) .GT.T)  GOTO  101 
TDELTA  =  T  -  TFTIME (KTFID) 

IF  (TDELTA. EQ. 0. )  THEN 
TFEXP  =  1. 

TFEXOT  =  l./TAUl (KTFID) 

GOTO  1 
END  IF 

TFEXP  =  EXP (-TDELTA/TAU1 (KTFID) ) 

C  FOLLOWING  CODE  ASSURES  CORRECT  LIMITING  VALUES  AS  TAU  ->  0. 

IF  (ABS (TAUl (KTFID) /TDELTA)  . LT .  l.E-4)  THEN 

TFEXOT  =  1.0/TDELTA 
ELSE 

TFEXOT  =  TFEXP/TAUl (KTFID) 

END  IF 

1  DRIVE  =  XIN 

OMEGA  =  DRIVE 

KTF  =  KTFTYP (KTFID) 

IF  (KTF.EQ.4)  OMEGA=  ( 1 . -TAU2 (KTF ID ) /TAUl (KTF ID ) ) * 

+  DRIVE 

IF(KTF.EQ.l  .OR.  KTF.EQ.4)  XTMP  = 

+  TFAUX (KTFID) *TFEXP+ ( 1 . -TFEXP ) *OMEGA 
IF  (KTF.EQ.l)  XOUT  =  XTMP 
IF  (KTF.EQ.2)  THEN 

XTMP  =  TFEXOT* (DRIVE-TFAUX2 (KTFID) ) +TFAUX (KTFID) *TFEXP 
XOUT  =  XTMP 
XTMP2=  DRIVE 
END  IF 

IF(KTF.EQ.3)  THEN 

TEMP  =  TFAUX2 (KTFID) -DRIVE 

XTMP=TFAUX (KTFID) +TAU1 (KTFID) *TEMP* ( 1 . -TFEXP ) +DRIVE*TDELTA 
XOUT  =  XTMP 

XTMP2  =  TEMP*TFEXP+DRIVE 
END  IF 

IF  (KTF.EQ'.4)  XOUT  =  DRIVE*  (TAU2  (KTFID) /TAUl  (KTFID)  ) 

+  +  XTMP 

'F  (KTF.EQ.5)  THEN 

KK  =  TAU3 (KTFID) 

II  =  TAUl (KTFID) 

BETA  =  TAU2 (KTFID) 

FORCE  =  DRIVE 

DELZ  =  TFAUX (KTFID) 

DELDZ  =  TFAUX2  (KTFID) 

TLAST  =  TFTIME (KTFID) 

CALL  SOLAG (II , BETA, KK, DELZ, DELDZ, FORCE, T, TLAST, DELT,DELTD) 
XTMP  =  DELT 

XTMP2  =  DELTD 

XOUT  =  XTMP 

END  IF 
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o  o  o  o 


100  TFTIME (KTFID)  =  T 

TFAUX(KTFID)  =  XTMP 

TFAUX2 (KTFID)  =  XTMP2 

XOUT  =  XOUT*GAIN (KTFID) 

C  ABOVE  SAVES  XD  FOR  2ND  ORDER  LAG/OSCILLATOR  AND  DX/DT  FOR 
C  FIRST  ORDER  LAG  WITH  DIFFERENTIATOR.  FOR  USE  AS  INITIAL 
C  CONDITIONS  ON  NEXT  INTEGRATION  STEP. 

RETURN 

101  WRITE(8,200)  T, TFTIME (KTFID) , KTFID 
STOP 

102  WRITE(8,201)  KTFID 
STOP 

103  WRITE (8,202)  KTF 

STOP 

200  FORMAT ( '  ***ERROR***  IN  TRNF .  TIME= ' , FIO . 5 , '  IS  GREATER  THAN 

+  'LAST  TIME='  ,F10. 5, '  FOR  CHANNEL '  , 1 3 , //) 

201  FORMAT ( '  ***ERROR***  IN  TRNF.  CHANNEL  NO.  KTFID=',I4, 

+  '  IS  OUT  OF  RANGE.'//) 

202  FORMAT ( '  ***ERROR***  IN  TRNF.  TYPE  CODE  NUMBER  KTFID=',I4, 

+  '  IS  OUT  OF  RANGE.'//) 

END 

SUBROUTINE  SOLAG (II , BETA , KK , DEL Z , DELDZ , FORCE , T , TLAST , DELT , 

+  DELTD) 

DIFFERENTIAL  EQUATION 

II*yDDOT  +  BETA*YDOT  +  KK*Y  =  FORCE 
NOTE:  THE  CALLING  ROUTINE  MUST  SAVE  TLAST, DELZ  AND  DELDZ 
FOR  THE  NEXT  CALL  TO  SOLAG 
REAL  NUM  ,LAMB1,  LAMB2,  LAMB,  II,  KK 
IF  (BETA.LT.O.  .OR.  KK . LE . 0 . )  THEN 
WRITE (8,1000)  BETA,KK 
STOP 
END  IF 

IF  (II.LT.O.)  THEN 
WRITE  (8,1001)  II 
STOP 
END  IF 

BOV2I  =  BETA/(2.*II) 

RDCND  =  B0V2I**2  -KK/II 
FOVK  =  FORCE/KK 
IF (RDCND  .GT.  0 . ) GO  TO  100 
IF(RDCND  .EQ.  0 . ) GO  TO  120 
C  OSCILLATORY  SOLUTION 

LAMB  =-BOV2l 
OMEGA  =  SQRT  (-RDCND) 

NUM  =  OMEGA* (DELZ-FOVK) 

DENOM  =  DELDZ-LAMB* (DELZ-FOVK) 

LAMB  =  -BOV2I 

IF(NUM  .EQ.  0.  .AND.  DENOM  .EQ.  0.)  GO  TO  20 
PHI  =  ATAN2 (NUM, DENOM)  -  OMEGA*TLAST 
GO  TO  30 
20  CONTINUE 


22 


PHI  =  0. 

30  CONTINUE 

ETEMP  =  EXP (-LAMB*TLAST) 

IF(DELZ  -  FOVK  .EQ.  0 . ) GO  TO  40 

A=(DELZ  -  FOVK) *ETEMP/SIN (OMEGA*TLAST+PHI ) 

GO  TO  50 
40  CONTINUE 

A=DELDZ*ETEMP/(LAMB*SIN (OMEGA*TLAST  +  PHI) 

1  +OMEGA*COS (OMEGA*TLAST  +  PHI)  ) 

50  CONTINUE 

FTEMP  =  EXP(LAMB*T) 

DELT=  A*FTEMP*SIN (OMEGA*T  +  PHI)  +  FOVK 
DELTD  =  LAMB*A*FTEMP*SIN (OMEGA*T+PHI )  + 

1  OMEGA*A*FTEMP*COS {OMEGA*T+PHI ) 

GO  TO  150 
100  CONTINUE 

C  EXPONENTIAL  SOLUTION 

LAMBl  =  -BOV2I+SQRT (RDCND) 

LAMB2  =  -BOV2I-SQRT (RDCND) 

A  =  (DELZ  -  DELDZ/LAMB2  -  FOVK) *EXP (-LAMBl*TLAST) / 

1  (1.  -  LAMB 1/LAMB 2) 

B  =  (DELZ  -  DELDZ/LAMBl  -  FOVK) *EXP (-LAMB2*TLAST) / 

1  (1.  -  LAMB 2/LAMB 1) 

DELT  =  A*EXP (LAMB1*T)  +  B*EXP (LAMB2*T)  +  FOVK 
DELTD  =  LAMB1*A*EXP (LAMB1*T) +LAMB2*B*EXP (LAMB2*T) 

GO  TO  150 
120  CONTINUE 

C  CRITICALLY  DAMPED  SOLUTION 

LAMB  =  -BOV2I 

A= ( (DELZ-FOVK) * ( 1+TLAST*LAMB) -TLAST*DELDZ ) *EXP (-LAMB*TLAST ) 
B= (DELDZ+LAMB* (FOVK-DELZ) ) *EXP (-LAMB*TLAST) 

DELT  =  (A+B*T) *EXP (LAMB*T) +FOVK 

DELTD  =  (LAMB*A+B+LAMB*B*T) *EXP (LAMB*T) 

150  CONTINUE 
TLAST  =  T 
DELZ  =  DELT 
DELDZ  =  DELTD 
PLTURN 

1000  -ORMAT(/,'  ***-ERR  IN  SOLAG-***  BETA= ' , E12 . 4 , ' , KK= ' , E12 . 4 , / , 
+  '  ONLY  POSITIVE  DEFINITE  VALUES  ALLOWED.’,/) 

1001  FORMAT(/,'  ***-ERR  IN  SOLAG-***  II=',E12.4,/ 

+  '  NEGATIVE  VALUES  NOT  ALLOWED.',/) 

END 
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SAMPLE  DATA  INPUT  CARDS 


K  K 

T 

T 

T 

G 

X 

T  T 

A 

A 

A 

A 

• 

F  F 

U 

U 

U 

I 

• 

I  T 

1 

2 

3 

N 

m 

D  Y 

• 

• 

• 

.  P 

•  • 

• 

• 

• 

• 

• 

• 

• 

• 

0101 

.02 

0. 

0. 

4.0 

0.0 

0202 

.2 

0. 

0. 

2.0 

0.0 

0303 

.02 

0. 

0. 

2.0 

0.0 

0404 

.02 

.01 

0.0 

1.0 

0.0 

0504 

.02 

0. 

0.0 

1.0 

0.0 

0605 

1. 

3. 

.25 

1.0 

10. 

0705 

1. 

1. 

.25 

1.0 

10. 

0805 

1. 

.5 

.25 

1.0 

10. 

0905 

1. 

.0 

.25 

1.0 

10. 

1005 

0. 

.5 

.25 

1.0 

10. 

1105 

1. 

.5 

.0 

1.0 

10. 

1205 

.001 

.0002 

20. 

1.0 

10. 

1305 

.0002 

0.0 

5000  . 

1.0 

10. 

D  A 

X  A 


0.  PITCH  ATTITUD  .  HOLD 

0.  PITCH  SYNTHETIC  DAMPING 

0.  YAW  ATTITUDE  HOLD 

0.  YAW  LOS  RATE  FILTER 

0.  YAW  LOS  RATE  FILTER 

0.  PURE  DAMPING 

0.  CRITICALLY  DAMPED 

0.  DAMPED  OSCILLATION 

0.  PURE  OSCILLATION 

0. 

0. 

0. 

0. 


0 

0 

0 

0 

0 

0 

0 

0 
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SAMPLE  OUTPUT 


TRANSFER  FUNCTIONS  DEFINITIONS: 


KTFID  KTFTYP  TAUl  TAU2 


TAU3  GAIN 


X 


DX 


1  1  .02000  .00000 

DESCRIPTION:  PITCH  ATTITUDE  HOLD 


.00000  4.00000  .00000 

GAIN/(TAU1*S+1) 


.00000 


2  2 
DESCRIPTION: 


.20000  .00000  .00000 

PITCH  SYNTHETIC  DAMPING 


2.00000  .00000 

GAIN*S/(TAU1*S+1) 


.00000 


3  3 

.02000 

. 00000 

DESCRIPTION: 

YAW 

ATTITUDE 

HOLD 

4  4 

.02000 

.01000 

DESCRIPTION: 

YAW 

LOS  RATE 

FILTER 

.00000  2.00000  .00000  .00000 

GAIN/(S* (TAU1*S+1) ) 


.00000  1.00000  .00000  .00000 

GAIN* (TAU2*S+1)/(TAU1*S+1) 


5  4  .02000  .00000  .00000 

DESCRIPTION:  YAW  LOS  RATE  FILTER 

***ERROR  IN  TRNF  10***  TAU  OUT  OF  RANGE. 

SUBSTITUTING  1ST  ORDER  LAG.  KTFTYP  ->  1. 
***THIS  TRANSFER  FUNCTION  WILL  BE  TREATED  AS: 

5  1  .02000  .00000  .00000 

DESCRIPTION:  YAW  LOS  RATE  FILTER 


1.00000  .00000  .00000 
GAIN* (TAU2*S+1) /(TAUl *S+1) 


1.00000  .00000  .00000 
GAIN/(TAU1*S+1) 


6  5 

DESCRIPTION: 

1.00000  3.00000 

PURE  DAMPING 

.25000 

1.00000  10.00000  .00000 
GAIN/(TAU1*S**2+TAU2*S+TAU3) 

7  5 

DESCRIPTION: 

1.00000 

CRITICALLY 

1.00000 

DAMPED 

.25000 

1.00000  10.00000  .00000 
GAIN/(TAU1*S**2+TAU2*S+TAU3) 

8  5 

DESCRIPTION: 

1.00000 
DAMPED  OSC 

.50000 

ILLATION 

.25000 

1.00000  10.00000  .00000 
GAIN/(TAU1*S**2+TAU2*S+TAU3) 

9  5 

DESCRIPTION: 

1.00000  .00000 
PURE  OSCILLATION 

.25000 

1.00000  10.00000  .00000 
GAIN/{TAU1*S**2+TAU2*S+TAU3) 
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10  5  .00000  .50000  .25000 

DESCRIPTION: 

***ERROR  IN  TRNF  10***  TAU  OUT  OF  RANGE. 

SUBSTITUTING  1ST  ORDER  LAG.  KTFTYP  ->  1. 
***THIS  TRANSFER  FUNCTION  WILL  BE  TREATED  AS: 
10  1  2.00000  .00000  .00000 

DESCRIPTION: 


1.00000  10.00000  .00000 

GAIN/(TAU1*S**2+TAU2*S+TAU3) 


4.00000  10.00000  .00000 

GAIN/(TAU1*S+1) 


11  5  1.00000  .50000  .00000  1.00000  10.00000  .00000 
DESCRIPTION:  GAIN/ (TAU1*S** 2+TAU2*S+TAU3 ) 
***ERROR  IN  TRNF  10***  TAU  OUT  OF  RANGE. 


SUBSTITUTING  1ST  ORDER  LAG  WITH  INTEGRATOR.  KTFTYP  ->  3. 
***THIS  TRANSFER  FUNCTION  WILL  BE  TREATED  AS: 


11  3 

2.00000 

.00000 

.00000 

2.00000  10.00000 

.00000 

DESCRIPTION: 

GAIN/(S* (TAU1*S+1) ) 

12  5  .00100 

DESCRIPTION: 


.00020  20.00000  1.00000  10.00000  .00000 

GAIN/(TAU1*S**2+TAU2*S+TAU3) 


13  5  .00020 

DESCRIPTION: 


.000005000.00000  1.00000  10.00000  .00000 

GAIN/{TAU1*S**2+TAU2*S+TAU3) 
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THE  INTEGRATION  TIME  STEP  DT  =  .005000  . 


LTFID 

TIME 

DRIVE 

X 

1 

.00000 

.00000 

.00000 

1 

1.00000 

-1.41615 

-5.53438 

1 

2.00000 

-1.65364 

-6.71729 

1 

3.00000 

-.03983 

-.20403 

1 

4.00000 

-1.14550 

-4.44207 

1 

5.00000 

-1.83907 

-7.42804 

LTFID 

TIME 

DRIVE 

X 

2 

.00000 

.00000 

.00000 

2 

1.00000 

-1.41615 

-3.67251 

2 

2.00000 

-1.65364 

1.68669 

2 

3.00000 

-.03983 

2.25946 

2 

4.00000 

-1.14550 

-3.56729 

2 

5.00000 

-1.83907 

.70957 

LTFID 

TIME 

DRIVE 

X 

3 

.00000 

.00000 

.00000 

3 

1.00000 

-1.41615 

-1.04245 

3 

2.00000 

-1.65364 

-4.69789 

3 

3.00000 

-.03983 

-6.27757 

3 

4.00000 

-1.14550 

-6.97196 

3 

5.00000 

-1.83907 

-10.47893 

LTFID 

TIME 

DRIVE 

X 

4 

.00000 

.00000 

.00000 

4 

1.00000 

-1.41615 

-1.39987 

4 

2.00000 

-1.65364 

-1.66648 

4 

3.00000 

-.03983 

-.04542 

4 

4.00000 

-1.14550 

-1.12801 

5.00000 

-1.83907 

-1.84804 

LTFID 

TIME 

DRIVE 

X 

5 

.00000 

.00000 

.00000 

5 

1.00000 

-1.41615 

-1.38360 

5 

2.00000 

-1.65364 

-1.67932 

5 

3.00000 

-.03983 

-.05101 

5 

4.00000 

-1.14550 

-1.11052 

5 

5.00000 

-1.83907 

-1.85701 
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LTFID 

TIME 

DRIVE 

X 

6 

.00000 

.00000 

10.00000 

6 

1.00000 

-1.41615 

9.35272 

6 

2.00000 

-1.65364 

8.09988 

6 

3.00000 

-.03983 

7.04001 

6 

4.00000 

-1.14550 

6.35853 

6 

5.00000 

-1.83907 

5.38481 

LTFID 

TIME 

DRIVE 

X 

7 

.00000 

.00000 

10.00000 

7 

1.00000 

-1.41615 

8.97774 

7 

2.00000 

-1.65364 

6.33926 

7 

3.00000 

-.03983 

3. 54-^  39 

7 

4.00000 

-1.14550 

1.78784 

7 

5.00000 

-1.83907 

.10284 

LTFID 

TIME 

DRIVE 

X 

8 

.00000 

.00000 

10.00000 

8 

1.00000 

-1.41615 

8.82387 

8 

2.00000 

-1 . 65364 

5.36682 

8 

3.00000 

-.03983 

1.15018 

8 

4.00000 

-1.14550 

-1.86486 

8 

5.00000 

-1.83907 

-4.23518 

LTFID 

TIME 

DRIVE 

X 

9 

.00000 

.00000 

10.00000 

9 

1.00000 

-1.41615 

8.62980 

9 

2.00000 

-1.65364 

3.87712 

9 

3.00000 

-.03983 

-3.25255 

9 

4.00000 

-1.14550 

-9.90176 

9 

5.00000 

-1.83907 

-15.20981 

LTFID 

TIME 

DRIVE 

X 

10 

.00000 

.00000 

40.00000 

10 

1.00000 

-1.41615 

23.29664 

10 

2.00000 

-1.65364 

11.23015 

10 

3.00000 

-.03983 

5.73197 

10 

4.00000 

-1.14550 

2.81734 

10 

5.00000 

-1.83907 

-1.12071 

LTFID 

TIME 

DRIVE 

X 

11 

.00000 

.00000 

20.00000 

11 

1.00000 

-1.41615 

19.86679 

11 

2.00000 

-1.65364 

18.71997 

11 

3.00000 

-.03983 

16.91363 

11 

4.00000 

-1.14550 

15.57969 
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11 

5.00000 

-1.83907 

13.85090 

LTFID 

TIME 

\ 

DRIVE 

X 

12 

.00000 

.00000 

10.00000 

12 

•  1.00000 

-1.41615 

-9.10857 

12 

2.00000 

-1.65364 

8.06504 

12 

3.00000 

-.03983 

-7.32883 

12 

4.00000 

-1.14550 

6.51484 

12 

5.00000 

-1.83907 

-5.97223 

LTFID 

TIME 

DRIVE 

X 

13 

.00000 

.00000 

10.00000 

13 

1.00000 

-1.41615 

1.54640 

13 

2.00000 

-1.65364 

-9.52188 

13 

3.00000 

-.03983 

-4.49206 

13 

4.00000 

-1.14550 

8.13177 

13 

5.00000 

-1.83907 

7.00721 
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PROGRAM  USAGE 
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The  program  consists  of  three  subroutines  (TRNF  10,  TRNF,  and  SOLAG)  and  a 
main  program,  TFMAIN.*  TRNF  10  is  the  initialization  routine.  The  data  records 
that  define  the  transfer  functions  are  read  by  this  routine.  (See  the  explana¬ 
tion  for  input  data  format  given  below.)  TRNF  performs  the  analytical  transfer 
function  simulations.  It  accesses  subroutine  SOLAG  for  treatment  of  the  second 
order  lag/harmonic  oscillator.  The  driving  main  program  TFMAIN  is  supplied  for 
demonstration  and  testing  only.  Use  it  as  a  model  for  interfacing  TRNF  10,  TRNF, 
and  SOLAG  with  your  computer  program.  The  initial  time  is  at  line  100,  and  the 
integration  time  step  is  at  line  200.  The  driving  term  is  defined  in  the  lines 
following  300. 


The  call  to  TRNF  contains  four  arguments : 

DRIVE  Driving  function 

XOUT  Output  of  the  transfer  function  channel 

T  Time 

J  Channel  I.D.  number  LTFID 

One  input  data  record  is  required  for  each  transfer  function  variable  to  be 
integrated.  These  records  have  the  following  form: 


Column 

Format 

Variable 

Definition 

• 

1 

12 

LTFID 

I.D  number  of  channel;  values  are  from 

1  to  50 

- 

3 

12 

LTFTYP 

Code  for  type  of  transfer  function; 
values  are  1  to  5 

10 

F9.0 

TAUl 

First  constant,  for  all  LTFTYP 

19 

F9.0 

TAU2 

Second  constant,  for  LTFTYP  =  4  or  5 

28 

F9.0 

TAU3 

Third  constant,  for  LTFTYP=5 

37 

F9.0 

XO 

Initial  condition  for  variable 

46 

F9.0 

DXO 

Initial  condition  for  derivative 
of  variable 

55 

A26 

AA 

Alphanumeric  label 

• 

LTFID  is  an 
transfer  function 

integer 

channel 

variable  from  1  to  50  that  is  associated  with  each 
(output  variable)  to  be  simulated.  For  example,  a  lag 

m 

*  A  listing  of  a 
appendix  A. 

program 

implementing 

these  techniques  and  a  sample  run  are  in 
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might  be  applied  to  both  a  fin  yaw  deflection  command  and  to  a  fin  pitch  deflec¬ 
tion  command.  Each  of  these  two  channels  is  defined  as  a  separate  data  input 
record  and  is  assigned  an  identifying  integer  LTFID.  If  the  user  chooses  an 
integer  that  was  previously  assigned  on  another  input  data  record,  execution 
terminates  with  an  explanatory  message. 

LTFTYP  is  an  integer  variable  in  the  range  1  to  5  that  defines  the  type  of 
transfer  function: 


LTFTYP  Description  and  FORTRAN  code 

1  First-order  lag 

1/(TAU1*S  +  1) 


First-order  lag  with  differentiator 
S/(TAU1*S  +  1) 


3  First-order  lag  with  integrator 
1/[S*(TAU1*S  +  1)] 

4  Combined  lead/lag 
(TAU2*S  +  1)/(TAU1*S  +  1) 

5  Second-order  lag/harmonic  oscillator 
1/(TAU1*S-^  +  TAU2*S  -i  TAU3) 


TAUl  is  the  first  constant  (see  LTFTYP  above).  All  five  transfer  function 
types  (LTFTYP  =  1  to  5)  require  this  quantity  to  be  positive. 

TAU2  is  the  second  constant.  Lead/lag  and  second-order  lag/oscillator 
transfer  functions  (LTFTYP  =  4,5)  require  a  positive  value. 

TAU3  is  the  third  constant.  The  second-order  lag/oscillator  (LTFTYP  =  5) 
requires  a  positive  value. 

XO  is  the  initial  condition  on  the  output  variable,  i.e.,  y^  =  y  at  t  =  t^. 

DXO  is  the  initial  condition  on  the  derivative  of  the  output  variable  v'  = 
dy/dt  at  t  =  t^.  '  ^ 

AA  is  an  alphanumeric  label  of  26  characters  that  identifies  the  nature  of 
the  transfer  function  when  an  echo  of  the  input  data  is  printed  out  by  TRNF  10. 
Examples  might  be  "PITCH  SYNTHETIC  DAMPING"  or  "YAW  ATTITUDE  HOLD." 

If  an  input  variable  is  out  of  range,  the  program  prints  an  error  message 
and  reformulates  the  transfer  function  into  an  equivalent  form.  For  example,  if 
an  attempt  is  made  to  run  an  oscillator  (LTFTYP  =  5)  with  TAUl  =  0,  subroutine 
TRNF  10  will  substitute  the  equivalent  simple  lag  (LTFTYP  =  1),  as  can  be  seen  in 
the  sample  output  for  KTFID  =10.  Two  other  examples  appear  in  the  sample  output 
in  the  sample  execution,  i.e.,  for  KTFID  =  5  and  KTFID  =  11. 
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SYMBOLS 


Algebraic  Expressions 


"  A 

Constant  of  integration 

B 

Constant  of  integration 

D 

Coefficient  of  the  Laplace  operator  s  in  the  transfer  function  of 
the  second-order  lag/oscillator 

I 

Coefficient  of  the  square  of  the  Laplace  operator  s  in  the  transfer 
function  of  the  second-order  lag/oscillator 

K 

Constant  term  in  the  transfer  function  of  the  second— order 
lag/ oscillator 

L 

Damping  coefficient  -D/21  for  the  second-order  lag/oscillator 

s 

Laplace  transform  operator 

t 

Time ,  independent  variable 

T 

Time  constant 

D 

Driving  term 

if 

y 

Output  of  transfer  function,  dependent  variable 

yo 

Initial  condition  on  y  at  t  =  t 

y' 

Initial  condition  on  dy/dt  at  t  =  t^ 

z 

Auxiliary  dependent  variable  for  first-order  lag  with  differentia¬ 
tion  or  for  first-order  lag  with  integration 

0) 

Frequency  of  oscillation  for  second-order  lag/oscillator 

<p 

Phase  angle  of  oscillation  for  second-order  lag/oscillator 

y 
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FORTRAN  Variable  Names 


* 


AA 

Alphanumeric  label  for  printout  (80  characters)  ' 

4 

DX 

Initial  condition  on  derivative  when  KTFTYP=5 

KTFID 

I.D.  number  of  channel  or  variable  (range  is  1  to  50) 

KTFTYP 

Integer  variable  used  to  indicate  type  of  transfer  function  (range 

1  to  50)  ® 

T 

KTFTYP  =  1  is  first-order  lag 

Transfer  function:  GAIN/(TAU1*S+1) 

KTFTYP  =  2  is  first-order  lag  with  differentiation 

Transfer  function:  GAIN*S/(TAU1*S+1) 

KTFTYP  =  3  is  first-order  lag  with  integration 

Transfer  function:  GAIN/ [S*(TAU1*S+1) ] 

KTFTYP  =  4  is  combined  lead/lag 

Transfer  function:  GAIN*(TAU2*S+1/(TAU1*S+1 ) 

KTFTYP  =  5  is  second-order  lag/harmonic  oscillator 

Transfer  function:  GAIN/ ( TAUl *S* *2  +  TAU2*S  +TAU3) 

Time 

TAUl 

First  time  constant 

TAU2 

Second  time  constant  (Tj^  for  KTFTYP=4  or  D  for  KTFTYP) 

TAU3 

Third  time  constant  (K  for  KTFTYP=5) 

X 

Initial  condition  on  variable 

XIN 

Driving  term 

XOUT 

Output  of  transfer  function 

r 
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